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Abstract 

The  research  objectives  are  to  develop  effective  procedures  for  computing  fluid  flow  at 
ultra  high  pressures  where  fluids  exhibit  very  different  thermodynamic  behavior  than  the 
perfect  gas  and  incompressible  fluid  models  that  are  commonly  used  in  CFD  simulations. 
Three  issues  to  be  addressed  include:  developing  RANS  algorithms  for  arbitrary  fluid 
applications;  developing  efficient  properties  evaluation  procedures  for  arbitrary  fluids; 
and  extending  hybrid  RANS-LES  algorithms  to  high  pressures.  All  three  of  these  issues 
have  been  demonstrated.  A  generalized  fluid  model  that  is  independent  of  the  equation  of 
state  has  been  developed  and  demonstrated  for  a  wide  variety  of  flows.  An  adaptive 
properties  evaluation  method  that  is  competitive  in  time  with  perfect  gas  calculations  but 
allows  highly  complex  EOS  routines  like  REFPROP,  tabular  properties  data,  or  simple 
algebraic  relations  (such  as  Peng-Robinson)  has  been  developed.  Results  for  hybrid 
RANS-LES  models  that  take  advantage  of  the  arbitrary  scaling  used  for  convergence  and 
accuracy  control  in  the  general  EOS  formulation  have  also  been  demonstrated.  The 
hybrid  method  is  efficient  at  lower  cell  aspect  ratios  below  100.  Additional  work  is 
needed  to  provide  the  desired  robustness  at  cell  aspect  ratios  of  1000  and  higher.  The 
general  equation  of  state  method  has  been  applied  to  a  number  of  supercritical  fluid 
applications  including  regenerative  heat  transfer  in  a  rocket  engine  chamber  and 
preliminary  solutions  for  combustion  problems  in  which  the  fuel  or  oxidizer  is 
supercritical  prior  to  burning.  In  addition  the  general  equation  methodology  has  been 
applied  to  Maxwell’s  equations  to  provide  a  hyperbolic  time-marching  method  that 
applies  in  either  wave-like  Maxwell  regimes,  or  in  diffusion-like  MHD  regimes. 

Introduction 

Numerous  Air  Force  T&E  applications  use  fluids  at  ultra  high  pressures  where  complex 
equations  of  state  are  required.  CFD  solutions  at  these  pressures  require  that  the 
equations  be  expressed  in  a  form  applicable  to  liquids,  vapors  and  supercritical  regions. 
The  numerical  algorithm  must  be  applicable  over  all  thermodynamic  regimes  and  across 
all  Mach  number  ranges.  The  representative  thermodynamics  at  these  conditions  are 
described  by  advanced  EOS  techniques  such  as  REFPROP1  which  has  currently  been 
updated  in  a  companion  effort.2  The  foundations  of  the  numerical  algorithms  were 
developed  previously  based  upon  an  appropriate  scaling  of  the  equations,3  but  have  been 
updated  in  the  present  work.  An  efficient  properties  evaluation  procedure  has  been  be 
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devised  to  ensure  effective  incorporation  in  CFD  codes.  Finally,  high  Reynolds  numbers 
applications  require  extensions  of  LES  turbulence  techniques  to  these  general  fluids. 

Equation  Formulation 

The  general  equations  of  motion  for  multi-component,  multi-phase  fluids  can  be 
expressed  in  the  divergence  form,  V*y  =  0,  where  V  is  a  four-dimensional  space-time 
operator  and  J  -  Qet  +  Eex  +  Fey  +  Gez .  The  temporal  flux  vector  includes  the  density, 

velocity  components,  pressure  and  stagnation  enthalpy,  Q  =  (p,  pu,  pv,  pwy  ph°  -  pf  , 
with  similar  terms  for  the  spatial  fluxes.  Numerical  solutions  do  not  allow  the  use  of 
density  as  an  independent  variable  in  thermodynamic  ranges  where  density  is  essentially 
constant  (liquids),  so  Q  cannot  be  used  as  the  primary  dependent  variable.  This  difficulty 
can  be  circumvented  by  introducing  a  pseudo-time  term,  r,  for  time-marching  and 

choosing  the  primitive  variables,  Qp  =  (p,w,v, w,T)7,  as  the  solution  variable  so  that  the 
equation  system  is  expressed  as,  rdQpdz  +  V  •  &  =  0.  The  pseudo-time  coefficient,  T, 

provides  dimensional  consistency  and  depends  on  the  equation  of  state  through  the 
property  derivatives  pp,  pj ,  hpy  and  hj  which  can  be  obtained  from  the  Gibbs 

function.2  Appropriate  scaling  of  the  equations  indicates  that  efficient  convergence  and 
uniform  solution  accuracy  over  all  thermodynamic  and  Mach  number  regimes  may  be 
achieved  by  replacing  the  physical  property,  ppy  by  an  artificial  property,  pp  that  is 

based  upon  the  Mach,  Reynolds,  Strouhal  and  Froude  numbers  of  the  problem. 

Properties  Evaluation  Procedure 

Having  set  the  equations  of  motion  and  the  algorithm,  efficient  properties  evaluation  can 
be  obtained  by  using  a  Cartesian  adaptive  table  look-up  procedure  that  enables  user- 
defined  accuracy  over  the  thermodynamic  domain  of  interest.  The  Cartesian  adaptive 
method  results  in  a  tree  structure  that  provides  search  procedures  in  highly  non-uniform 
tables  that  are  competitive  with  equally  spaced  tables  while  also  providing  the  high 
accuracies  in  both  the  thermodynamic  functions  and  their  derivatives  that  are  required  for 
CFD  calculations.  Representative  results  for  CO2  in  Fig.  1  show  the  grid  for  a  0.1% 
accurate  solution.  The  evaluation  efficiency  of  this  method  in  REFPROP  tables  shows  a 
speed  increase  between  two  to  three  orders  of  magnitude  and  CFD  computations  based  on 
the  arbitrary  fluid  require  essentially  the  same  CPU  time  as  those  based  on  a  perfect  gas. 
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Figure  1.  Adaptive  interpolation  maps  for  C02.  10*2  <  p  <  103  MPa ,  500  <  T  <  1600  K. 
Left:  1%  accuracy,  22,000  points;  Right:  0.1%  accuracy,  225,000  points. 

Three  levels  of  interpolation  accuracy  have  been  considered:  C°,  C1  and  C2.  As  is 
discussed  in  more  detail  below,  at  all  three  levels  of  accuracy  the  final  square  domains  are 
broken  into  triangles  over  which  bi-variate  interpolation  is  accomplished.  Cell 
boundaries  across  which  a  change  in  refinement  level  result  in  property  discontinuities 
that  can  adversely  impact  CFD  convergence  and  accuracy  unless  care  is  taken  to  generate 
triangles  that  are  consistent  with  the  grid  refinement.  For  the  C°  interpolation,  the  two 
thermodynamic  functions  and  their  derivatives  are  all  taken  as  independent  linear 
functions.  This  implies  the  derivatives  are  not  consistent  with  the  functions,  but  ensures 
derivative  continuity  across  interpolation  cell  boundaries.  The  C1  interpolation  uses  fifth- 
order  interpolation  to  give  second-order  continuity  in  density  and  enthalpy  at  all 
triangular  boundaries.  The  enthalpy  and  its  derivatives  are  then  evaluated  consistently  in 
each  triangle  and  are  continuous  at  cell  boundaries.  The  density  and  its  derivatives  are 
similarly  consistent  and  continuous.  The  C2  interpolation  uses  ninth-order  interpolation 
to  give  second-order  continuity  in  the  Gibbs  function  so  that  the  enthalpy,  density  and  all 
four  derivatives  are  mutually  consistent  and  continuous.  Results  show  that  all  three 
methods  provide  effective  CFD  solutions. 

Thermodynamic  Consistency  and  Function  Continuity  in  Properties  Reconstruction 

For  a  fluid  computation  based  upon  primitive  variables,  the  six  thermodynamic 
properties,  p ,  pp,pT,  h,hp,  and  hj ,  and  the  two  transport  properties,  p  and  k,  must  be 
evaluated  as  a  function  of  p  and  T  at  every  cell  and  every  iteration.  There  are  several 
potential  ways  to  employ  Cartesian  adaptive  methods  to  store  and  recover  the  data.  The 
most  straightforward  approach  is  to  store  all  eight  properties  at  each  node  of  the  Cartesian 
grid  and  reconstruct  all  eight  properties  independently.  While  this  procedure  can  provide 
continuous  functions  across  cell  boundaries,  it  produces  inconsistent  results  because  it 
ignores  the  interdependency  of  the  density  and  the  enthalpy  and  the  four  thermodynamic 
derivatives.  To  incorporate  this  characteristic  in  CFD  solutions,  the  consistency  of  the 
property  evaluations  must  be  considered. 


Thermodynamic  consistency  in  property  evaluations  has  been  discussed  by  Swesty,5  who 
used  the  Helmholtz  free  energy  as  a  basis  function.  His  results,  which  were  based  on 
much  coarser  meshes  (and  correspondingly  lower  accuracies)  than  the  ones  contemplated 
here,  show  that  failing  to  include  proper  consistency  conditions  in  thermodynamic 
reconstruction  can  lead  to  significant  errors  in  flowfield  solutions.  Our  experience  with 
fine  grid  interpolation  normally  keeps  the  inconsistencies  to  acceptable  levels  so  that  we 
can  use  either  consistent  or  inconsistent  properties  evaluation.  Swesty  used  higher-order 
reconstruction  to  maintain  second-order  consistency  throughout  a  piece-wise  continuous 
domain  of  quadrilaterals.  In  analogous  fashion,  we  consider  three  different  levels  of 
reconstruction  involving  C°,  C\  and  C2  continuity  across  cell  boundaries  whose  results 
range  from  partially  to  fully  consistent.  In  place  of  the  quadrilateral  sub-domains  used  by 
Swesty,  we  reconstruct  the  property  functions  over  triangles  in  the  mapped  domain.  With 
the  global  structure  maintained  through  the  quadrilateral  grid,  properties  are  retrieved  by 
first  locating  the  triangle  in  which  their  given  p  and  T conditions  lie  and  then  using  stored 
values  at  the  vertices  to  interpolate  to  the  desired  degree  of  accuracy  during  a  numerical 
simulation. 

In  the  C°  method,  values  for  all  thermodynamic  and  transport  properties  are  stored  at  the 
vertices  of  each  square  in  the  mapped  plane.  The  squares  are  then  subdivided  into 
triangles  and  the  values  of  each  function  at  the  three  vertices  are  used  to  construct 
independent  bi-linear  reconstruction  functions  for  each  property.  The  bi-linear 
reconstruction  used  in  the  C°  method  therefore  ensures  function  continuity  along  all  faces 
of  adjacent  cells  of  the  same  size  for  all  eight  fluid  properties  but  does  not  provide 
internal  consistency  inside  the  triangles.  For  example,  both  density  and  its  first 
derivatives  vary  linearly  over  each  triangular  region.  Further,  the  enthalpy  and  density 
have  no  interdependence  on  each  other.  Nevertheless,  by  taking  the  six  properties  stored 
at  the  vertices  from  a  consistent  thermodynamic  database,  the  internal  inconsistencies 
remain  small  on  a  fine  grid. 

The  alternative  of  using  linear  reconstruction  over  all  triangles  for  the  density  and 
enthalpy  and  differentiating  these  linear  functions  to  obtain  the  property  derivatives  gives 
improved  consistency,  but  introduces  derivative  discontinuities  at  the  cell  faces. 
Experience  shows  that  these  derivative  discontinuities  are  much  more  detrimental  to  CFD 
solutions  than  the  inconsistencies. 

To  incorporate  property  consistency  in  the  reconstructed  results  we  use  bivariate 
polynomials  of  the  form22  23, 

i(x.y)=YJY^q,.,x'yi  (1) 

(=0  y=0 

to  represent  the  local  solution.  The  number  of  undetermined  constants  in  a  polynomial  of 
degree  n  is  (n  +  lXn  +  2)/2 .  One  advantage  of  the  Cartesian  subdivision  is  that  these 
polynomials  are  especially  suitable  for  triangles.  The  theory  and  implementation  of 
polynomials  on  triangular  meshes  is  well  documented  (see,  for  example,  Akimo6  and 
Preusser,7)  and  both  C1  and  C2  reconstruction  methods  based  upon  the  polynomial 
expansion  of  Eq.  1  have  been  considered.  The  C1  reconstruction  provides  consistency 
between  density  and  its  partial  derivatives  and  enthalpy  and  its  partial  derivatives,  but 
does  not  provide  consistency  between  density  and  enthalpy.  We  refer  to  the  C1 


reconstruction  as  partially  consistent.  The  C2  reconstruction  provides  fully  consistent 
results.  Details  of  both  methods  are  given  below. 

The  C1  method  treats  the  density  and  enthalpy  as  unrelated  functions  that  are  each 
handled  in  analogous  fashion.  Using  the  density  as  an  example,  the  goal  of  C1 
reconstruction  is  to  reconstruct  the  density  function  in  such  a  manner  that  the  values  of 
p ,  pp  and  pT  are  consistent  within  any  given  cell  and  that  all  three  quantities  are 
continuous  across  cell  boundaries  (i.e.,  the  density  is  C1  continuous  while  the  derivatives 
are  C°  continuous).  This  level  of  consistency  and  continuity  can  be  achieved  by  using  a 
bi-quintic  (fifth-order)  polynomial  of  the  form  given  in  Eq.  1 .  The  construction  of  a  bi- 
quintic  polynomial  requires  the  determination  of  21  coefficients.  Storing  the  values  of 
the  density  and  its  derivatives  up  to  second  order  at  each  of  the  three  vertices  provides  18 
conditions.  The  remaining  three  come  from  the  requirement  that  the  derivatives  in  the 
direction  normal  to  each  edge  be  continuous.  This  can  be  done  either  by  prescribing  the 
normal  derivatives  or  by  reducing  the  order  of  the  polynomials.  The  later  choice  is  the 
known  as  the  ‘condensation  of  parameters’  method.6  Since  the  density  and  enthalpy  are 
treated  independently,  the  enthalpy  and  its  first  two  derivatives  must  also  be  stored  at  all 
vertices.  Similar  consistency  and  continuity  are  then  ensured  for  the  enthalpy  and  its  first 
derivatives.  The  enthalpy  and  density  are,  however,  unrelated  on  a  given  region  so  that 
consistency  between  these  two  functions  is  not  guaranteed.  Nevertheless,  if  both  density 
and  enthalpy  are  taken  from  a  database  such  as  REFPROP  that  is  itself  consistent,  the 
degree  of  inconsistency  in  the  Cl  approximation  is  small. 

The  C2  method  interpolates  all  properties  and  their  derivatives  by  reconstructing  the 
Gibbs  function  as  a  ninth-order  polynomial  which  requires  the  determination  of  55 
coefficients.  Of  these,  45  are  determined  by  storing  the  Gibbs  function  and  its  derivatives 
up  to  fourth  order  on  the  three  vertices.  Requiring  (^-continuity  on  each  of  the  edges 
provides  another  nine  conditions.  The  final  condition  can  be  specified  at  will  without 
violating  the  smoothness  constraints,  while  still  keeping  the  fit  to  the  nodal  data. 
Additional  details  of  the  formulation  are  given  in  the  literature.6  In  the  C2  method,  only 
the  Gibbs  function  is  reconstructed  and  all  six  thermodynamic  properties  are  computed 
from  this  single  function  resulting  in  properties  that  are  fully  consistent  with  each  other 
and  the  underlying  Gibbs  function.  All  six  thermodynamic  properties  are  based  upon 
higher  order  variations  of  the  Gibbs  function  across  each  triangle  with  first-  and  second- 
derivative  continuity  across  all  adjoining  faces.  This  results  in  consistency  and  continuity 
among  all  six  properties.  The  resulting  density  and  enthalpy  reconstructions  are 
consistent  and  C1  continuous,  while  their  derivatives  are  consistent  and  Cf  continuous. 

For  all  three  reconstruction  levels,  the  coefficients  of  the  piecewise  function  are  stored  for 
each  triangle  once  they  have  been  calculated  and  need  not  be  updated  in  CFD 
calculations.  In  the  following  section  we  discuss  appropriate  triangulation  methods  for 
effectively  implementing  the  C°,  C\  and  C2  methods  on  the  adaptive  Cartesian  grid. 
Triangulation  on  Unequally  Sized  Squares.  In  equally  spaced  regions  of  a  Cartesian 
mesh  where  all  neighbors  are  the  same  size,  the  squares  can  be  split  into  triangles  by 
dividing  along  either  diagonal  and  any  of  the  three  reconstruction  methods  described 
above  can  be  applied  on  the  resulting  triangles.  The  situation  is  somewhat  more  complex 
when  two  neighboring  cells  are  unequal.  For  example,  the  line  /  shown  in  Fig.  2  which 
crosses  face  AF  between  cells  of  two  different  refinement  levels.  The  simple  diagonal 


triangulation  indicated  in  Fig.  2  results  in  a  property  discontinuity  at  this  interface  for 
either  the  (?,  C1  or  C2  reconstruction  procedures  because  the  properties  on  the  left  side  of 
the  interface  (point,  Pi)  are  determined  by  values  of  the  function  at  vertices  A  and  F, 
while  those  on  the  right  side  (point  PR)  are  determined  by  vertices  A  and  D. 


Figure  2.  Details  of  diagonal  triangulation  on  a  square  mesh  at  a  change  in  cell  size. 

This  difficulty  is  demonstrated  in  Fig.  3  which  compares  the  reconstructions  along  line  / 
of  the  density  (left)  and  its  pressure  derivative,  pp  (right),  based  upon  the  triangulation 
shown  in  Fig.  2  with  the  exact  solution  as  functions  of  temperature.  In  addition,  the 
reconstruction  based  upon  an  alternative  triangulation  discussed  below  is  also  presented. 
The  discontinuity  at  the  interface  for  the  simple  diagonal  triangulation  (Fig.  4)  is  readily 
visible.  The  reconstructions  in  Fig.  3  were  obtained  by  using  the  C°  method,  but 
reconstructions  based  on  C1  or  C2  exhibit  similar  discontinuities.  Further,  the 
discontinuity  exists  regardless  of  which  diagonal  is  used  to  bisect  the  squares. 


Figure  3.  Reconstructed  values  of  p  and  pp  for  triangulations  in  Fig.  2  (Long  dashed  line) 

and  Fig.  6  (short-dashed  and  dotted  lines)  along  line  /  in  Figs.  2  and  6.  Exact  solution  given  by 
dotted  points. 


To  demonstrate  the  impact  that  discontinuities  of  this  nature  have  on  CFD  calculations, 
we  show  the  results  of  one-dimensional  CFD  solution  on  Fig.  4  and  the  corresponding 
convergence  plots  on  Fig.  5.  The  axial  variation  of  the  temperature  in  Fig.  4  indicates 
that  the  discontinuity  in  properties  at  a  cell  size  change  produces  an  unphysical  wiggle  in 
the  solution.  The  corresponding  convergence  curve  on  Fig.  5  also  indicates  that 
convergence  stalls  when  it  reaches  a  level  consistent  with  the  magnitude  of  the 
discontinuity  in  the  reconstructed  functions.  The  remaining  convergence  curves  are 
discussed  below. 


Figure  4.  Temperature  variation  in  1-D  computational  solution  with  reconstruction  based 
on  diagonally  divided  squares  (Fig.  2)  illustrating  wiggles  created  in  solution  by  property 
discontinuities  at  cell  faces. 

Modified  Triangulation  for  Ensuring  Continuity.  The  difficulties  arising  from 
unequally  sized  squares  can  be  rectified  by  triangulating  the  Cartesian  mesh  in  a  manner 
that  ensures  the  reconstructions  on  both  sides  of  the  interface  are  based  on  the  same  data. 
One  acceptable  triangulation  pattern  for  a  cell  with  two  refined  and  two  unrefined 
neighbors  is  demonstrated  in  Fig.  6.  As  shown  by  the  subdivision  in  the  figure,  the 
properties  on  the  left  and  right  sides  of  the  boundary  at  point  P  are  each  based  upon 
similar  values  and  result  in  the  specified  degree  of  function  continuity.  The 
reconstructions  of  the  properties  along  line  l  across  the  boundary  are  now  both 
determined  by  the  values  at  A  and  F.  The  continuity  of  the  piecewise  function  across  the 
interface  is  thus  guaranteed  as  exhibited  by  the  (?,  C 1  and  C2  reconstruction  curves  in 
Fig.  5.  The  C°  reconstruction  results  in  a  small  local  error,  but  this  has  no  adverse  impact 
on  the  convergence  rate  (Fig.  7)  or  the  solution  (not  shown).  The  reconstructed  functions 
for  the  C 1  and  C2  reconstruction  also  provide  excellent  convergence  as  Fig.  7  shows. 
Finally,  Fig.  7  shows  the  convergence  rate  for  a  solution  for  which  the  properties  were 
taken  directly  from  REFPROP  without  any  reconstruction  method.  This  ‘exact’  property 
evaluation  method  resulted  in  exactly  the  same  convergence  as  any  of  the  three 


reconstruction  levels  (although  it  took  much  more  CPU  time).  Consequently,  we  see  that 
a  proper  triangulation  results  in  both  efficient  convergence  and  accurate  solutions. 


1 


-10 

-12 


- Square  (Ditcentfnous) 

T  riangl*(C°-Co  ntinou*) 

-  Quintic  (C1 -C  orttinout) 

-  Nonic  (C? -C ontinou*) 

•  Exact  Solution  (Rafprop) 


_ _ _ _ _ _ _ I _ _ _ _ _ _ _ I - - - - - - - 1 - - - . - . - 1 - - - i - - - 1 

0  200  400  600  800  1000 

Iterations 


Figure  5.  Exemplary  plot  showing  effect  of  triangulation  on  convergence  of  a  specific  CFD 
calculation.  Red:  Triangulation  taken  from  Fig.  4;  Green,  blue  and  violet:  Triangulation 
taken  from  Fig.  6;  Triangles:  properties  evaluated  directly  from  REFPROP. 


Figure  6.  Triangulation  method  chosen  to  ensure  continuous  reconstruction  in  a  region 
between  two  cell  sizes. 

Similar  triangulations  can  be  found  for  nearly  any  type  of  mesh  topology,  but  to  limit  the 
number  of  patterns  to  a  countable  level,  we  prohibit  cell  level  changes  greater  than  one  in 
neighboring  cells  during  the  initial  database  setup.  Thus  the  ratio  of  face  length  in  two 
adjacent  cells  will  be  either  one  or  two.  This  limitation,  which  is  also  commonly  used  in 
tree-based  adaptive  grid  methods,16  17  reduces  the  number  of  patterns  to  six:  cells  with 
four  undivided  faces,  cells  with  one  divided  and  three  undivided  faces;  cells  with  two 


adjacent  faces  divided;  cells  with  two  opposite  faces  divided;  cells  with  three  faces 
divided;  and  cells  with  all  four  faces  divided.  These  six  patterns  are  shown  in  Fig.  7 
along  with  appropriate  triangulations.  Because  only  six  patterns  must  be  recognized,  the 
logic  for  the  triangulation  is  straightforward  and  the  overhead  for  locating  the  triangles 
within  a  given  square  is  small  so  that  the  search  advantage  of  the  Cartesian  mesh  is  not 
lost.  The  reconstruction  process  applies  equally  to  any  triangle  whether  it  is  a  portion  of  a 
regular  region  of  the  mesh  or  one  with  unequal  sides.  Consequently,  the  modified 
triangulation  adds  no  complexity  to  the  reconstruction  procedure  and  allows  the  overall 
structure  of  the  Cartesian  grid  to  be  retained. 


Figure  7.  Six  possible  cell  refinement  patterns  indicating  acceptable  triangulation  patterns  for 
each. 


Storage  and  Timing  Comparisons  Next  we  present  some  timing  comparisons  for  the 
adaptive  reconstruction  procedure  and  some  representative  reconstruction  results.  We 
consider  properties  evaluations  for  three  fluids,  CO2,  H2O  and  air.  Statistically 
significant  information  for  timing  is  obtained  by  considering  several  different  zones  in  the 
p-T  domain  for  each  fluid.  Some  of  the  reconstruction  regions  include  the  discontinuity 
across  the  liquid-vapor  line  while  some  are  restricted  to  continuous  property  regions.  For 
these  two-phase  regions,  we  also  compare  thermodynamic  procedures  that  treat  the  two 
phases  as  a  single  fluid  with  procedures  that  treat  them  as  two  distinct  fluids.  Following 
the  definition  of  these  various  zones,  we  show  the  tree-structure  grid  for  reconstruction  to 
both  1%  and  0.1%  accuracy.  Finally  we  compare  the  time  required  for  the  table  look-up 
procedure  and  the  complete  REFPROP  solution. 

Figure  8  shows  the  density  of  C02  as  a  function  of  pressure  and  temperature  as  obtained 
from  REFPROP.  The  liquid-vapor  line  is  shown  along  with  the  global  density  variation. 
Four  different  rectangular  zones  are  used  for  reconstruction  databases.  Zone  1  comprises 
the  temperature  and  pressure  ranges,  500 K  <  T  <  1 600/f  ,  0.01  MPa  <  p  <  \  000MPa  and 
lies  entirely  within  the  vapor  region.  The  grid  structures  for  this  zone  are  shown  on  Fig.  9 
for  accuracies  of  1%  and  0.1%.  The  grid  color  is  keyed  to  the  magnitude  of  the  density. 
The  1%  accuracy  case  requires  a  total  of  ten  levels  of  refinement  and  results  in  22,000 
reconstruction  points.  The  0.1%  case  requires  1 1  refinement  levels  and  225,000  points. 
Similar  results  are  obtained  for  Zones  2  and  3.  The  table  fit  and  reconstruction  are  done 
on  the  basis  of  the  logarithm  of  the  pressure.  The  size  of  the  reconstruction  maps  is 
summarized  in  Table  I. 
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Figure  8.  Density  of  C02  as  a  function  of  pressure  and  temperature  (from  REFPROP). 
Reconstruction  databases  have  been  obtained  for  four  zones:  Zones  1  and  2  do  not  include  the 
vanor-liauid  line:  Zones  3  and  4  overlao  the  discontinuity. 


Figure  9.  Adaptive  reconstruction  maps  for  Zone  1  of  C02  (see  Fig.  10).  Left:  1%  accuracy, 
22,000  points;  Right:  0.1%  accuracy,  225,000  points.  Timing  evaluations  made  by  performing 


Zones  3  and  4  of  the  CO2  map  cross  the  liquid  vapor  line  where  the  density  is 
discontinuous.  Maps  with  accuracies  of  1%  and  0.1%  are  shown  on  Fig.  10.  The 
presence  of  the  discontinuity  implies  that  the  refinement  criterion  can  never  be  satisfied  at 
the  liquid-vapor  line.  We  have  arbitrarily  terminated  the  refinement  process  for  this  case 


at  15  levels.  This  implies  that  the  errors  will  be  less  than  the  stated  values  at  all  locations 
except  those  immediately  adjacent  to  the  liquid-vapor  line. 


Figure  10.  Adaptive  reconstruction  maps  for  Zone  4  (Fig.  10)  of  C02.  Left:  1%  accuracy,  72,703 
points;  Right:  0.1%  accuracy,  108,384  points.  Domain  location:  220/f  <  7  <  500A"  ; 


An  alternative  reconstruction  map  for  the  CO2,  Zone  4  domain,  is  given  in  Fig.  1 1  but 
with  the  vapor  and  liquid  regions  broken  into  two  separate  curve  fits.  Because  it  is  no 
longer  necessary  to  attempt  to  interpolate  across  the  discontinuity,  the  number  of 
reconstruction  cells  for  a  given  accuracy  is  reduced  dramatically.  At  the  1%  accuracy 
level,  the  number  of  cells  reduces  from  22,000  for  the  continuous  fit  across  the 
discontinuity  to  5336  points  in  the  vapor  and  225  points  in  the  liquid.  Thus  a  total  of 
5561  points  (and  only  ten  levels  of  refinement)  give  the  same  (or  slightly  improved) 
accuracy  in  the  ‘two-phase’  description  as  in  the  ‘single-phase’  description.  For  the  0.1% 
accuracy  map,  the  number  of  reconstruction  cells  is  reduced  from  108,000  to  52,296  in 
the  vapor  plus  1 102  points  in  the  liquid  region. 


Figure  11.  Adaptive  interpolation  maps  for  Zone  4  of  C02  (see  Fig.  10).  Left:  1%  accuracy,  225 
points  in  liquid  zone,  5336  points  in  vapor  zone;  Right:  0.1%  accuracy,  1102  points  in  liquid  zone. 


Similar  properties  databases  have  been  generated  for  water  and  air.  The  density  contour 
maps  for  these  two  fluids  are  shown  in  Figs.  12  and  13  along  with  two  zones  for  which 
properties  data  have  been  generated  in  a  manner  similar  to  that  discussed  above  for  CO2. 
For  brevity,  the  Cartesian  maps  for  these  cases  are  not  given  here,  but  the  resulting  table 
sizes  a  d  the  reconstruction  times  are  included  in  Table  1  discussed  below. 
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Figure  12.  Density  of  H20  as  a  function  of  pressure  and 
temperature  (taken  from  REFPROP).  Interpolation  databases 


Figure  13.  Density  of  Air  as  a  function  of  pressure  and  temperature  (taken 
from  REFPROP).  Interpolation  databases  have  been  obtained  for  the  two 
zones  indicated. 


CPU  Comparisons  with  C°  Reconstruction  Method.  An  assessment  of  the  property 
evaluation  times  required  for  the  adaptive  reconstruction  procedure  as  compared  to  that 
for  the  original  REFPROP  routine  is  given  in  the  bar  chart  on  Fig.  14.  Timing 
comparisons  are  shown  for  four  zones  of  carbon  dioxide,  two  zones  of  water  and  two 
zones  of  air.  The  evaluation  times  are  plotted  on  a  logarithmic  scale  with  three  bars  given 
for  each  fluid  zone.  The  bars  represent,  respectively,  the  time  for  reconstruction 
accuracies  of  1%  and  0.1%  along  with  the  time  for  the  REFPROP  evaluations.  All 
timings  are  based  on  making  100,000  property  evaluations  at  equally  spaced  points  on 
diagonal  lines  across  the  respective  fluid  zones  (as  indicated  by  the  dashed  lines  on  Fig.  9 
for  CO2).  The  reconstruction  timings  are  for  the  C°  method.  The  added  cost  for  the  C1 
and  C2  methods  is  indicated  later. 


Figure  14.  Comparison  of  CPU  time  of  adaptive  Cartesian  property  reconstruction 
and  REFPROP  calculation. 

The  timing  evaluations  in  Fig.  14  immediately  show  the  advantage  of  the  reconstruction 
procedure.  In  all  cases,  reconstruction  is  at  least  two  orders  of  magnitude  faster  than  the 
complete  REFPROP  routines.  The  evaluation  times  for  reconstruction  are  essentially  the 
same  for  all  fluid  types,  all  fluid  regions,  and  for  the  two  accuracy  levels.  This 
insensitivity  of  evaluation  time  to  accuracy  level  (table  size)  is  a  clear  indication  of  the 
effectiveness  of  the  tree  structure  in  the  reconstruction  table.  By  contrast,  the  evaluation 
times  for  the  REFPROP  routines  vary  considerably  with  fluid  type  and  fluid  region 
because  the  complexity  of  the  underlying  equations  changes.  Air,  which  is  treated  as  a 
mixture,  is  considerably  more  expensive  to  evaluate  than  the  pure  fluids.  Because  of  the 
timing  variability  in  REFPROP,  the  smallest  time  ratio  occurs  in  Zone  1  for  H20  where 
reconstruction  is  85  times  faster  while  the  maximum  advantage  occurs  for  air  where 


reconstruction  is  as  much  as  2500  times  faster.  The  savings  achieved  by  the 
reconstruction  procedure  would  clearly  be  smaller  for  equations  of  state  based  upon 
simpler  algebraic  equations,  but  even  for  the  simplest  equations,  the  procedure  remains 
competitive. 
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Figure  15.  Storage  Comparison  of  different  reconstruction  methods.  C02, 
1%  accuracy.  Domain  location:  500 K  <  T  <  1000/(  ;  \MPa  <  p  <  100 MPa  . 


The  timing  results  in  Fig.  14  are  also  tabulated  in  Table  I  along  with  the  number  of 
refinement  levels  and  the  total  number  of  cells  in  each  zone.  The  times  required  for 
making  the  100,000  property  evaluations  are  also  given  for  an  equally  spaced  table  to 
indicate  the  penalty  incurred  by  the  tree-structure.  The  cost  of  a  tree-based  reconstruction 
is  approximately  12%  larger  than  that  for  the  equally  spaced  database.  The  overhead 
arises  mainly  from  the  coordinate  transformation  and  local  triangle  selection  within  a 
square.  A  numerical  value  of  the  ratio  of  the  adaptive  reconstruction  time  to  the 
REFPROP  evaluation  time  is  also  included.  The  columns  in  Table  I  are  defined  in  Table 

n. 

One  caution  in  interpreting  these  comparisons  is  that  REFPROP  uses  the  Helmholtz 
function  which  expresses  fluid  properties  as  functions  of  temperature  and  density  as  the 
fundamental  thermodynamic  variable.  Property  evaluations  in  REFPROP  based  upon 
pressure  and  temperature  will  therefore  incur  some  penalty.  To  demonstrate  that  this  is 
not  the  major  reason  for  the  advantage  of  the  adaptive  method,  the  CO2  properties 
evaluation  in  Zone  2  was  chosen  as  a  representative  comparison.  The  time  required  for 
100,000  property  evaluations  at  equally  spaced  points  based  upon  pressure  and 
temperature  was  76.5  s  as  noted  in  Table  I.  Similar  evaluations,  based  upon  density  and 
temperature  required  45.1  s,  a  savings  of  42  percent  over  the  pressure  and  temperature 
pair.  Therefore  for  CFD  codes  that  are  based  upon  density  and  temperature,  the  CPU 
advantages  of  the  adaptive  Cartesian  method  would  be  reduced  by  approximately  a  factor 
of  two,  still  providing  a  major  savings.  For  codes  based  upon  other  variables  pairs,  the 
speed-ups  in  Table  I  remain  appropriate. 


Table  I.  Storage  and  Timing  Comparisons  for  Properties  Evaluation 


Fluid 


co2 


h2o 


Air 


Zone 

Max 

No. 

No. 

Time(s) 

Time(s) 

Time(s) 

Time 

Err 

Levels 

Points 

Cartes 

E.E.M* 

REFPROP 

Ratio 

1% 

10 

22,246 

0.438 

0.375 

77.686 

177.4 

1 

0.1% 

11 

225,121 

0.438 

0.375 

77.297 

176.5 

1% 

17 

25,990 

0.422 

0.375 

76.532 

181.4 

2 

0.1% 

20 

25,7841 

0.422 

0.375 

76.250 

180.7 

1% 

15* 

15,9325 

0.422 

0.375 

73.022 

173.1 

3 

0.1% 

15* 

80,1151 

0.438 

0.375 

73.250 

167.2 

1% 

15* 

72,703 

0.422 

0.375 

87.329 

207.1 

4 

0.1% 

15* 

108,384 

0.422 

0.375 

87.563 

207.5 

4 

1% 

5 

86 

(Liquid) 

0.1% 

7 

1102 

4 

1% 

10 

5336 

(Vapor) 

0.1% 

12 

52,296 

1 

1% 

10 

28,391 

0.438 

0.375 

36.234 

82.7 

1 

0.1% 

12 

225,121 

0.438 

0.375 

36.343 

83.0 

1% 

9 

871 

0.406 

0.375 

186.107 

458.4 

2 

0.1% 

11 

8049 

0.422 

0.375 

186.796 

442.6 

1% 

11 

63.740 

0.422 

0.375 

202.672 

480.3 

1 

0.1% 

12 

79,3657 

0.453 

0.375 

202.203 

446.4 

1% 

8 

459 

0.391 

0.375 

989.203 

2530.0 

2 

0.1% 

10 

4184 

0.406 

0.375 

989.532 

2437.3 

The  lowest  level  allowed  (for  cases  where  discontinuity  presents) 
Equivalent  Equal-sized  Mesh  (A  uniform  mesh  refined  to  the  deepest  level) 


Table  II  Column  Information  in  Table 


Column 

Heading 

Information 

1 

Fluid 

Fluid  Type 

2 

Zone 

Zone 

3 

Max  Err 

Maximum  error  of  curve  fit 

4 

No.  Levels 

Number  of  refinement  levels  needed  (or  allowed  in  cases  with  property 
discontinuity) 

5 

No.  Points 

Number  of  points  in  the  adaptive  database 

6 

Cartes  Time 

Time  required  for  100,000  property  evaluations  using  adaptive  database 

7 

Time(s)  E.E.M. 

Time  required  for  100,000  property  evaluations  using  equally  spaced 
database. 

8 

REFPROP  Time 

Time  required  for  100,000  property  evaluations  using  REFPROP 

9 

Time  Ratio 

Ratio  of  REFPROP  to  adaptive  database  property  evaluation  time 

Table  Size  Comparison  of  C1  and  C2  Reconstruction  Methods.  The  above  results 
have  shown  the  CPU  time  advantages  of  the  adaptive  reconstruction  method.  In  the 
present  subsection,  we  compare  the  table  sizes  required  for  the  Cfy  C 1  and  C2 
reconstruction  procedures.  The  comparisons  are  obtained  for  the  rectangular 
temperature-pressure  region,  500/f  <  T  <  1000/f  ,  \MPa  <  p  <  100 MPa,  in  CO2  using  an 
error  tolerance  of  1%.  .  In  the  lower  pressure  portion  of  this  region,  CO2  behaves  as  a 
perfect  gas,  while  in  the  upper  regimes  it  is  strongly  supercritical  (see  Fig.  8).  The  sizes 
of  the  storage  tables  for  the  three  interpolation  levels  are  given  in  Fig.  15.  In  computing 
these  tables,  the  Gibbs  function  and  its  first  two  derivatives  were  obtained  from  analytical 
expressions  within  REFPROP  but  the  high  order  derivatives  (greater  than  two)  were 
obtained  by  numerical  differentiation.  The  number  of  cells  required  for  the  C° 
reconstruction  method  was  1290,  while  the  C1  method  decreased  the  table  size  to  only  23 
cells.  The  C2  reconstruction,  however,  required  75  cells  to  achieve  the  same  accuracy, 
suggesting  a  sensitivity  in  the  C2  reconstruction  procedure  to  the  high  order  derivatives. 
This  phenomenon  may  arise  from  two  issues.  First,  the  Gibbs  function,  rather  than  the 
density  and  enthalpy  themselves,  is  fitted  in  the  C2  reconstruction  and  the  density  and 
enthalpy  are  obtained  by  differentiating  the  Gibbs  function.  Second,  the  C2 
reconstruction  uses  higher  order  polynomials  that  are  more  prone  to  oscillate  and  this 
tendency  may  require  increased  resolution  to  realize  a  given  accuracy. 

To  demonstrate  that  this  trend  in  the  storage  size  is  not  related  to  the  method  in  which  the 
REFPROP  evaluation  are  made,  we  have  added  in  Fig.  15  companion  reconstruction 
results  for  a  perfect  gas  using  the  same  pressure-temperature  region  in  CO2.  Here,  all 
partial  derivatives  were  computed  analytically  and  the  specific  heats  were  taken  as 


constant.  The  perfect  gas  tables  for  Cfy  C1  and  C2  reconstruction  have  664,  14  and  42 
mesh  points  respectively,  again  indicating  a  significant  savings  in  moving  from  Cf] 
reconstruction  to  its  C1  counterpart  with  a  modest  increase  for  C2  reconstruction. 

Sample  CFD  Calculations 

In  the  present  section  we  present  some  representative  computational  results  to 
demonstrate  the  method.  We  begin  with  a  one-dimensional  computation  that  is  used  to 
verify  the  method  and  to  provide  timings.  We  then  present  a  two-dimensional  example 
that  demonstrates  some  of  the  thermodynamic  characteristics  that  are  enabled  by  the 
present  formulation.  These  computations  are  based  on  the  results  of  an  in-house  CFD 
code.  The  code  uses  a  second-order  accurate,  approximate  Riemann  solver  in  space  with 
a  two  equation  turbulence  model.  The  focus  in  these  examples  is  on  demonstrating  the 
efficiency  and  practicality  of  real  fluid  computations. 

One-Dimensional  Flow  Computation.  As  an  initial  example  of  the  timing  realized  with 
the  adaptive  reconstruction  method,  we  use  a  simple  one-dimensional  example  of 
subsonic  flow  through  a  convergent-divergent  nozzle.  The  nozzle  is  symmetric  about  the 
throat  with  the  inlet  and  exit  areas  1.25  times  the  throat  area.  The  working  fluid  is  chosen 
as  water  with  an  inlet  total  pressure  of  60  A/Pa  and  stagnation  temperature  of  750  K .  The 
back  pressure  at  the  outlet  is  50  MPa.  These  parameters  place  the  solution  in  Zone  2  of 
the  H20  property  map  on  Fig.  14.  A  total  of  200  cells  are  used  in  the  calculation.  The 
convergence  rates  for  this  case  using  the  C°,  Cl  and  C2  methods  were  originally  shown  in 
Fig.  7  where  we  noted  that  the  convergence  with  all  three  of  these  reconstruction  methods 
was  identical  to  the  convergence  obtained  by  coupling  the  REFPROP  routines  directly 
into  the  CFD  code,  so  long  as  the  triangulation  depicted  in  Fig.  8  is  used. 
Correspondingly  we  showed  that  convergence  on  the  triangulation  of  Fig.  2  stalled. 
Although  the  number  of  iterations  is  essentially  identical  for  the  different  property 
evaluation  methods,  their  CPU  costs  are  significantly  different. 

Figure  16  shows  the  cumulative  CPU  time  required  for  the  computation  as  a  function  of 
iteration  number  for  1000  iterations.  The  (f  reconstruction  method  results  in  the  fastest 
execution.  Comparison  of  the  results  for  the  inconsistent  and  the  consistent  triangulation 
methods  indicates  that  the  subdivided  triangular  mesh  increases  the  CPU  time  by  about 
5%  as  compared  to  the  diagonally  divided  triangles,  indicating  a  minor  overhead  for 
locating  a  particular  triangle  inside  each  square.  The  costs  of  the  C1  and  C2  methods  are 
nearly  equal  and  are  approximately  twice  that  of  the  C°  method.  The  C2  method  evaluates 
all  properties  from  a  single  ninth-order  polynomial  from  while  the  C 1  method  evaluates 
properties  from  separate  fifth-order  polynomials  for  density  and  enthalpy.  The 
calculation  based  on  the  exact  property  evaluation  from  REFPROP  is  approximately  160 
times  slower  than  the  (f  method  and  80  times  slower  than  the  C1  and  C2  methods,  again 
in  keeping  with  the  timing  of  the  properties  evaluation  themselves. 

The  CPU  time  for  a  CFD  calculation  can  be  broken  into  two  parts,  the  property 
evaluation  time  and  the  equation  solution  time.  The  relative  advantages  of  faster  property 
evaluation  times  clearly  will  decrease  as  the  complexity  of  the  equation  solution 
increases,  and  in  particular  as  we  move  from  one  to  three  dimensions.  The  nominal  cost 


of  the  property  evaluation  by  REFPROP  in  this  region  of  the  H20  map  is  about  450  times 
that  of  the  C°  integration  method  (see  Table  I).  For  this  one-dimensional  case,  the  CFD 
calculation  time  ratio  is  160.  suggesting  that  the  equation  solution  time  is  roughly  80 
percent  longer  than  the  Cf  interpolation  time  in  the  CFD  calculation.  Experiments  on  a 
PC  show  that  the  one-dimensional  solver  costs  about  6\  .25 {is /(cell  x  iteration) ,  and  our 
in-house  2D  and  3D  unstructured  codes  cost  about  100.2  {is/ (cell  x  iteration)  and  200.9 
{is / (cell  x  iteration)  respectively.  Assuming  the  property  evaluation  time  does  not 
change  when  going  from  ID  to  2D  and  3D  calculations,  this  indicates  that  the  2D 
calculation  with  the  C°  interpolation  method  will  be  115  times  faster  than  the  direct 
REFPROP  evaluation  calculation  and  the  3D  calculation  will  be  66  times  faster. 
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Figure  16.  CPU  time  comparison  for  different  interpolation  methods 


The  CPU  time  for  a  CFD  calculation  can  be  broken  into  two  parts:  the  property 
evaluation  time  and  the  equation  solution  time.  The  relative  advantages  of  faster  property 
evaluation  times  clearly  will  decrease  as  the  complexity  of  the  equation  solution 
increases,  and  in  particular  as  we  move  from  one  to  three  dimensions.  The  nominal  cost 
of  the  property  evaluation  by  REFPROP  in  this  region  of  the  H20  map  is  about  450  times 
that  of  the  Cf  integration  method  (see  Table  I).  For  this  one-dimensional  case,  the  CFD 
calculation  time  ratio  is  160,  suggesting  that  the  ratio  of  equation  solution  time  to 
property  evaluation  time  in  a  one-dimensional  calculation  is  approximately  1.8. 
Numerical  computations  show  that  the  one-dimensional  solver  requires  about 
61 .25 {is /(cell  x  iteration) .  Corresponding  costs  for  two-  and  three-dimensional  solutions 
are  100.2  and  200.9  {is/ (cell  x  iteration)  respectively.  Assuming  the  property  evaluation 
time  does  not  change  when  going  from  one  to  three  dimensions,  this  suggests  that  a  two- 
dimensional  calculation  with  Cf  interpolation  method  will  be  115  times  faster  than  the 
direct  REFPROP  evaluation  calculation  and  the  three-dimensional  calculation  will  be  66 
times  faster.  Solutions  based  upon  Cy  and  C2  reconstruction  will  be  about  half  this 
amount. 


Two-Dimensional  Real-Fluid  Example.  The  geometrical  configuration  for  the  2-D 
computation  is  again  a  C-D  nozzle  patterned  after  one  designed  for  hypersonic  flow 
testing  at  ultra  high  pressures.28  The  working  fluid  is  air  with  upstream  stagnation 
conditions  of  1700  MPa  and  750  K.  The  density  of  air  at  these  conditions  is 
approximately  equal  to  that  of  water  and  the  resulting  acceleration  through  the  choked 
throat  shows  dramatic  real-fluid  characteristics.  The  convergent  section  of  the  nozzle  is 
very  strongly  converging  and  the  inflow  starts  at  low  speeds. 

We  limit  our  results  here  to  a  Mollier  chart  for  air  showing  the  physical  regime  of  interest 
and  comparisons  between  solutions  of  the  real-fluid  and  corresponding  perfect  gas 
computations  at  the  same  pressure  and  temperature.  The  Mollier  diagram  is  shown  on  the 
left  half  of  Fig.  17  along  with  the  h-s  domain  covered  by  the  flow  in  the  nozzle. 
Specifically,  the  enthalpy  and  entropy  values  in  all  cells  in  the  computational  domain  are 
plotted  on  the  h-s  diagram  superimposed  upon  the  thermal  database  from  REFPROP. 
The  solution  forms  a  roughly  triangular  region.  The  left  edge  of  this  triangular  region 
corresponds  to  the  thermodynamic  path  of  the  fluid  on  the  centerline.  As  can  be  seen,  the 
fluid  on  the  centerline  undergoes  an  approximately  constant  entropy  process.  The  upper 
edge  of  the  triangle  corresponds  to  the  fluid  adjacent  to  the  wall  and  shows  its  entropy 
continues  to  increase  while  the  enthalpy  decreases  slightly  in  agreement  with  classical 
results  for  the  adiabatic  recovery  temperature  on  a  nozzle  wall.  The  horizontal  constant 
pressure  lines  to  the  right  side  of  the  Mollier  diagram  correspond  to  the  perfect  gas  region 
where  enthalpy  is  independent  of  pressure.  The  nearly  vertical  constant  pressure  lines  on 
the  left  side  indicate  a  strong  real-fluid  effect  of  pressure  on  enthalpy.  These  lead  to  a 
very  different  temperature  pattern  in  the  nozzle  as  the  contours  in  the  right  half  of  Fig.  17 
show. 
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Figure  17.  2D  calculation  of  ultra-high  pressure  flow  through  hypersonic  nozzle.  Left:  Mollier  diagram  and  h- 
s  domain  of  solution;  Right:  Temperature  along  the  nozzle  for  calculation  with  real  fluid  properties  and  perfect 

The  temperature  distribution  in  the  nozzle  for  a  perfect  gas  calculation  is  compared  with 
that  for  the  real-fluid  calculation  on  the  right  half  of  Fig.  17.  Here,  the  temperature  in 


every  computational  cell  has  been  plotted  at  its  appropriate  axial  location  in  the  nozzle. 
The  left  side  of  the  plot  shows  the  temperature  for  the  perfect  gas  solution  with  the  axial 
coordinate  running  from  left  to  right.  The  right  side  of  the  plot  shows  the  results  for  the 
real-fluid  case  with  the  axial  coordinate  running  from  right  to  left.  There  is  clearly  a 
major  difference  between  the  two  solutions.  In  both  solutions,  the  upper  bound  on  the 
temperature  at  any  axial  location  corresponds  to  the  wall,  while  the  minimum  temperature 
represents  the  value  in  the  free  stream.  For  the  perfect  gas,  the  temperature  on  the  wall 
(the  upper  bound)  decreases  slightly  from  its  upstream  value  of  750  K ,  while  for  the  real 
fluid,  the  wall  temperature  rises  rapidly  throughout,  reaching  over  1400  K  at  the  throat  (x 
=  0)  and  continuing  to  above  1800  K.  Contrasting  the  increased  wall  temperature  with 
the  nearly  constant  stagnation  enthalpy  for  the  real  fluid  calculation  in  the  Mollier 
diagram  on  the  left  plot  of  Fig.  17  demonstrates  that  the  wall  enthalpy  for  the  real  fluid 
behaves  like  that  for  the  perfect  gas  in  that  it  remains  essentially  constant  throughout  the 
nozzle.  The  wall  temperature  for  the  real-fluid  solution,  however,  increases  rapidly 
because  of  the  strong,  real-fluid,  pressure  dependence  of  enthalpy.  The  minimum 
temperature  at  any  axial  location  is  also  lower  for  the  real  fluid  than  for  the  perfect  gas, 
again  indicating  how  much  the  real-fluid  thermodynamics  differ  from  the  perfect  gas 
assumption.  As  a  final  point,  we  note  that  both  the  convergence  rates  and  the  CPU  times 
for  the  perfect  gas  and  real  fluid  computations  were  nearly  identical,  thus  verifying  that 
the  reconstruction  procedure  provides  an  effective  thermodynamic  interface  for  real-fluid 
computations. 

Hybrid  RANS/LES  Modeling  in  High  Pressure  Fluids 

Contemporary  turbulence  modeling  efforts  emphasize  large  eddy  simulations  (LES)  and, 
more  recently,  hybrid  RANS/LES  models,  but  efforts  in  these  areas  have  been  largely 
limited  to  incompressible  fluids  and  perfect  gases.  Application  of  these  methods  to 
arbitrary  fluids  requires  careful  adaptation  to  ensure  that  the  thermodynamics,  the 
numerics  and  fluid  dynamics  are  properly  coordinated  to  give  accurate  and  efficient 
results.  Our  focus  here  is  on  developing  methods  that  can  be  applied  to  practical  T&E 
applications  to  take  into  account  flow  fields  that  contain  both  traditional  high  Reynolds 
number  boundary  layers  and  regions  with  large  scale  unsteadiness  that  is  not  accurately 
resolved  by  unsteady  RANS  methods.  For  this  reason,  our  focus  is  on  hybrid  RANS/LES 
models  that  are  compatible  with  traditional  second-order  CFD  algorithms. 

To  adapt  the  generalized  fluid  dynamic  formulation  outlined  above  to  hybrid  RANS/LES 
modeling,  it  is  necessary  to  reduce  the  level  of  artificial  dissipation  in  the  algorithm.  The 
numerical  dissipation  in  characteristics-based,  upwind  methods  is  determined  by  the 
eigenvalues  of  the  problem  and  the  time-marching  derivative.  The  artificial  pseudo-time 
coefficient  matrix  offers  flexibility  for  controlling  this  dissipation.  The  flux  formulation 
for  a  general  upwind  finite  volume  approach  can  be  interpreted  as  the  average  of  the 
fluxes  on  either  side  of  the  cell  interfaces  augmented  by  an  artificial  dissipation  term, 


Clearly,  the  magnitude  of  the  artificial 


dissipation  is  proportional  to  the  coefficient  matrix,  r.  The  key  requirements  on  this 
term  are  that  it  provide  sufficient  dissipation  to  prevent  solution  instability  and  that  it 


scale  properly  at  various  flow  limits  including  low  speeds,  low  Reynolds  numbers  and 
high  Strouhal  numbers.  To  reduce  the  artificial  dissipation,  we  introduce  a  modified 
dissipation  matrix  that  improves  accuracy  while  maintaining  convergence,  so  that  the 


numerical  flux  becomes,  E  =  ^(EL  +  ER)-^r\D\(QpR  -  QpL)  where  D  is  a  diagonal 


matrix  in  which  one  entry  corresponds  to  the  largest  pseudo-acoustic  value  and  the 
remainder  are  the  particle  speed,  while  a  is  a  parameter  to  be  adjusted.  Analytical  and 
numerical  experiments  suggest  that  the  scheme  is  unstable  for  a>  1  and  becomes  stiff  as 
a  — »  0 .  The  value  a  =  0.5  represents  a  compromise  between  stability  and  accuracy  and 
is  used  in  the  calculations  below. 


The  modified  artificial  dissipation  is  first  verified  by  comparing  with  the  analytical 
solution  for  an  infinite  array  of  counter-rotating  vortices  whose  amplitudes  decay  with 
time.  The  results  on  the  left  of  Fig.  18  for  inviscid  flow  and  Re  =  100  were  computed 
with  CFL  =  0.1.  The  standard  (‘non-preconditioned’)  calculation  shows  that  both  cases 
decay  rapidly  demonstrating  the  accepted  observation  that  unmodified  upwind  schemes 
are  unacceptable  for  DES  applications.  When  ‘characteristic  dissipation’  is  used,  the 
inviscid  case  decays  by  about  1%  in  two  non-dimensional  time  units,  while  the  Re  =  100 
case  decays  by  4.7%.  Although  this  is  a  major  improvement,  the  decay  rates  are  still  too 
fast.  With  the  modified  dissipation,  the  inviscid  case  decays  by  approximately  0.04%  in 
two  non-dimensional  time  units  while  the  Re  =  100  case  decays  approximately  3.7%, 
slightly  slower  than  the  exact  solution.  Accurate  calculations  require  that  both  the  CFL 
and  the  VNN  numbers  be  controlled  and  this  small  error  becomes  larger  as  Re  is  reduced. 
To  demonstrate  accuracy  at  lower  Re  numbers  we  show  on  the  right  side  of  Fig.  18 
results  for  a  CFL  of  0.0036  for  Reynolds  numbers  of  infinity,  100,  10  and  1.  At  Re  =  1, 
CFL  =0.0036  corresponds  to  VNN  =  0.36,  while  at  CFL  =  0.1  and  Re  =  100  the  von 
Neumann  number  is  0.1.  The  results  at  this  viscous-controlled  time  step  indicate  the 
modified  dissipation  continues  to  provide  accuracy  at  all  Reynolds  numbers.  Overall,  the 
results  of  the  Taylor  vortex  problem  suggest  the  modified  dissipation  provides  a 
reasonable  solution  for  large-scale  motions. 


Figure  18.  Rate  of  decay  of  kinetic  energy  at  different  Reynolds  numbers  for  Taylor 
vortices  problem.  Left:  CFLU  =  0.1;  Right:  CFLU  =  0.0036 .  Comparison  with  analytical 
solution. 

The  scaling  constant,  CDES ,  in  the  DES  implementation  was  calibrated  by  comparing 
results  against  the  decay  of  homogeneous  turbulence  downstream  of  a  grid  as  reported  by 
Comte-Bellot  and  Corrsin.  Figure  19  shows  the  energy  spectra  for  calculations  using  the 
characteristic  and  modified  dissipation  terms  along  with  the  experimental  results  and  a 

line  showing  the  £(**)-  slope.  The  value  of  the  DES  coefficient,  CDES ,  in  these 
cases  is  0.78.  It  can  be  clearly  seen  that  the  modified  dissipation  significantly  improves 
the  simulation.  When  CDES  is  calibrated  against  experiment  (the  right  side  of  Fig.  19) 
the  turbulence  decay  is  commensurate  with  the  experiment.  In  Fig.  19  CDES  of  0.5 
appears  to  give  the  best  agreement  with  experiment  and  is  adopted  in  the  solutions  given 
below. 

Representative  DES  solutions  are  compared  with  unsteady  RANS  calculations  in  Fig.  420 
for  a  shear  layer  representing  the  mixing  between  air  and  nitrogen  at  high  pressures.  The 
results  indicate  substantially  better  resolution  for  the  modified  scheme.  Efforts  are 
present  are  aimed  at  extending  these  results  to  larger  grid  aspect  ratio  cases  (above  100) 
to  enable  extensions  to  higher  Reynolds  numbers. 


Figure  19.  Energy  spectra  in  decaying  homogeneous  turbulence  for  different 
preconditioning  methods  (left)  and  for  different  CDES  values  (right). 


Figure  20.  Comparison  of  RANS  (left)  and  DES  (right)  calculations  for  shear  layer  mixing 
in  high  pressure  flow. 
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